﻿namespace LxCFD
{
    public class PossionEquationBase : EquationBase
    {
        public Value1Base<double> Gamma;

        protected void setBound_flux0(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ALUO[i1] = ALUN[i1] = 0;
                b_FaceCommon[i1] = 0;
            }
        }
        protected void setBound_flux(IndexNumber indexNumber, ValueBase<double> flux, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            if (flux is Value1Base<double> flux1)
            {
                for (; i1 < endp1; i1++)
                {
                    ref FvFace face = ref Faces[i1];
                    ALUO[i1] = ALUN[i1] = 0;
                    double fluxValue = flux1.GetFaceValue(i1) * face.Area;
                    b_FaceCommon[i1] = fluxValue;

                    Vector3d v = face.Normal - face.gT * face.OToN;
                    double svv = Vector3d.DotProduct(face.Normal, v) / Cells[face.OSideCell].Volumn;
                    b[face.NSideCell] = fluxValue / Gamma.GetFaceValue(i1)
                        - Vector3d.DotProduct(Grad[face.OSideCell], v) + svv * x[face.NSideCell];
                    ADiag[face.NSideCell] = svv + face.gT;
                }
            }
            else if (flux is Value2Base<double> flux2)
            {
                for (; i1 < endp1; i1++)
                {
                    ref FvFace face = ref Faces[i1];
                    ALUO[i1] = 0;
                    flux2.GetFaceValue(out double sa, out double sb, i1);
                    sa *= face.Area;
                    sb *= face.Area;
                    b_FaceCommon[i1] = sb + sa * (x[face.NSideCell] - x[face.OSideCell]);
                    ALUN[i1] = sa;

                    Vector3d v = face.Normal - face.gT * face.OToN;
                    double svv = Vector3d.DotProduct(face.Normal, v) / Cells[face.OSideCell].Volumn;
                    double gamma = Gamma.GetFaceValue(i1);
                    b[face.NSideCell] = sb / gamma
                        - Vector3d.DotProduct(Grad[face.OSideCell], v) + svv * x[face.NSideCell];
                    ADiag[face.NSideCell] = svv + face.gT - sa / gamma;
                }
            }
        }
        protected void updateBound_flux(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                x[face.NSideCell] = (b[face.NSideCell] + face.gT * x[face.OSideCell]) / ADiag[face.NSideCell];
            }
        }
        protected void updateBound_fixhFixT(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                x[face.NSideCell] = (b[face.NSideCell] - x[face.NSideCell] * x[face.OSideCell]) / ADiag[face.NSideCell];
            }
        }
        protected void setBound_fixhFixT(IndexNumber indexNumber, Value1Base<double> fixh, Value1Base<double> fixT, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double hS = fixh.GetFaceValue(i1) * face.Area;
                double hS_div_gamma = hS / Gamma.GetFaceValue(i1);
                double PhiInf = fixT.GetFaceValue(i1);
                double FluxCb = hS * face.gDiff / (hS_div_gamma + face.gDiff);
                ALUN[i1] = -FluxCb;//会被加到ADiag
                ALUO[i1] = 0;
                b_FaceCommon[i1] = FluxCb * PhiInf + hS * face.CalNonOrthogonalityCorrection_bound(Grad, x) / (hS_div_gamma + face.gDiff);

                double xn_xo = x[face.NSideCell] - x[face.OSideCell];
                Vector3d v = face.Normal - face.gT * face.OToN;
                double svv = Vector3d.DotProduct(face.Normal, v) / Cells[face.OSideCell].Volumn;
                b[face.NSideCell] = hS_div_gamma * PhiInf
                    - (Vector3d.DotProduct(Grad[face.OSideCell], v) - svv * xn_xo);
                double atemp = svv + face.gT;
                ADiag[face.NSideCell] = atemp + hS_div_gamma;
                x[face.NSideCell] = -atemp;
            }
        }
        protected void setAb_step1_inner(int taskIndex)
        {
            InnerIndex.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                ALUO[i1] = ALUN[i1] = -gamma * face.gDiff;
                b_FaceCommon[i1] = gamma * face.CalNonOrthogonalityCorrection_inner(Grad, x);
            }
        }
        protected void setAb_step1_bound(IndexNumber indexNumber, int taskIndex)
        {
            indexNumber.FaceIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvFace face = ref Faces[i1];
                double gamma = Gamma.GetFaceValue(i1);
                ALUO[i1] = ALUN[i1] = -gamma * face.gDiff;
                b_FaceCommon[i1] = gamma * face.CalNonOrthogonalityCorrection_bound(Grad, x);
            }
        }
        protected void setAb_step2(Value2Base<double> source2,int taskIndex)
        {
            InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
            for (; i1 < endp1; i1++)
            {
                ref FvCell cell = ref Cells[i1];
                source2.GetCellValue(out double S_Ad, out double S_b, i1);
                ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO) - S_Ad * cell.Volumn;
                b[i1] = cell.CalIfOwnerAddElseMunus(b_FaceCommon) + S_b * cell.Volumn;
            }
        }
        protected void setAb_step2(Value1Base<double> source1,int taskIndex)
        {
            if (source1.IsEqual0())
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    ref FvCell cell = ref Cells[i1];
                    ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO);
                    b[i1] = cell.CalIfOwnerAddElseMunus(b_FaceCommon);
                }
            }
            else
            {
                InnerIndex.CellIDs.GetSeg(out int i1, out int endp1, taskIndex);
                for (; i1 < endp1; i1++)
                {
                    ref FvCell cell = ref Cells[i1];
                    ADiag[i1] = -cell.CalAllAdd(ALUN, ALUO);
                    b[i1] = cell.CalIfOwnerAddElseMunus(b_FaceCommon) + source1.GetCellValue(i1) * cell.Volumn;
                }
            }
        }
    
    }
}
